/*------------------<--- Start of Description -->--------------------\
| Returns Schoenfeld residuals and scaled Schoenfeld residuals from |
| the Cox model. The Schoenfeld residuals(r) are defined only at |
| event times. There is one residual for each covariate in the Cox |
| model. The scaled residuals(Bt) are defined as: |
| Bt= B + D*cov(B)*r', B=cox model beta, |
| D=total #events and |
| r=Schoenfeld residuals. |
| Bt is an estimate of the hazard function at follow-up time t. |
| Therefore a plot of Bt vs time can be used to assess whether the |
| proportional hazards assumption is valid. As time is frequently |
| skewed, plots of Bt vs rank time or '1-Kaplan-Meier' for the entire|
| dataset(which is a monotone function of time) may be preferred. |
| The correlation of Bt with time provides a test of the |
| proportional hazards assumption. The Breslow method is used for |
| tied event times. |
|--------------------<--- End of Description -->---------------------|
|--------------------------------------------------------------------|
|--------------<--- Start of Files or Arguments Needed -->-----------|
| Arguments: |
| - Required: |
| time=time variable for survival |
| event=event variable for survival, 1=event, 0=censored |
| xvars=list of covariates for Cox model |
| - Optional: |
| data= name of analysis dataset. Default is last dataset |
| created. |
| strata=stratification variable(one only) for stratifed Cox |
| models. |
| outsch=name of output data set containing Schoenfeld |
| residuals. One obs per each event time. The variables |
| containing the residuals have the same name as the |
| covariates (xvars). The data set also includes the time|
| variable and the strata variable. |
| Default data set name is schr. |
| outbt =name of output data set containing the scaled Schoenfeld|
| residuals(Bt). One obs per each event time. The |
| variables containing the scaled residuals have the same |
| name as the covariates (xvars). The dataset also |
| includes the time variable, the rank of the time(rtime) |
| and a variable(probevt) which is equal to '1 - overall |
| Kaplan-Meier' at the given time. |
| Default data set name is schbt. |
| plot=t,r,k,n. Indicates that SAS/Graph plots of Bt vs |
| time(t), rank of time(r) or '1 - overall Kaplan-Meier' |
| (k) are to be done. Default is r. For no plots use |
| n. The name of the graphics catalog is gschbt. |
| vref= indicator to control plotting of a vertical reference |
| line at y=0. Values are yes(default) and no. |
| points=yes,no. Whether to plot the actual data points. |
| Default is yes. |
| df= degrees of freedom for smoothing process Possible |
| values are 3 - 7. Default is 4. |
| pvars= variables to plot. Default is all xzvars. |
| alpha= confidence coefficient for plotting standard |
| error bars. Default is .05. A value of 0 means |
| do not plot SE bars. |
| rug= indicator to control plotting of rug of x values. |
| Values are yes and no(default). |
| Output: The macro prints the PHREG output used to fit the Cox |
| model, correlation coefficients of Bt vs time and |
| SAS/Graph plots of Bt vs time. |
|---------------<--- End of Files or Arguments Needed -->------------|
|--------------------------------------------------------------------|
|----------------<--- Start of Example and Usage -->-----------------|
| Usage: %schoen(data=_last_,time=,event=,xvars=,vref=yes, |
| strata=,outsch=schr,outbt=schbt,plot=r,points=yes, |
| df=4,pvars=,alpha=.05,rug=no); |
| References: Schoenfeld, D. (1982). Partial residuals for the |
| proportional hazards regression model. Biometrika 69, 239-41. |
| Grambsch PM, Therneau TM (1993). Proportional hazards tests |
| and diagnostics based on weighted residuals. University of |
| Minnesota, Division of Biostatistics. Research Report 93-002. |
\-------------------<--- End of Example and Usage -->---------------*/
%macro schoen(data=_last_,time=,event=,xvars=,vref=yes,
strata=,outsch=schr,outbt=schbt,plot=r,points=yes,
df=4,pvars=,alpha=.05,rug=no);
/*--------------------------------------------\
| Author: E. Bergstralh and T. Therneau; |
| Created: May 6, 1993; |
| Purpose: Returns Schoenfeld residuals and |
| scaled Schoenfeld residuals from |
| the Cox model; |
\--------------------------------------------*/
footnote"schoen macro: event=&event time=&time strata=&strata";
footnote2"Xvars= &xvars";
%let bad=0;
%if &time= %then %do;
%put ERROR: NO TIME VARIABLE GIVEN.;
%let bad=1;
%end;
%if &event= %then %do;
%put ERROR: NO EVENT VARIABLE GIVEN.;
%let bad=1;
%end;
%if &xvars= %then %do;
%put ERROR: NO XVARS GIVEN.;
%let bad=1;
%end;
%if &df<3 | &df>7 %then %do;
%put ERROR: DF MUST BE BETWEEN 3 AND 7.;
%let bad=1;
%end;
%if %upcase(&points)^=YES & %upcase(&points)^=NO %then %do;
%put ERROR: POINTS MUST BE YES OR NO.;
%let bad=1;
%end;
%if %upcase(&rug)^=YES & %upcase(&rug)^=NO %then %do;
%put ERROR: RUG MUST BE YES OR NO.;
%let bad=1;
%end;
%if %upcase(&vref)^=YES & %upcase(&vref)^=NO %then %do;
%put ERROR: VREF MUST BE YES OR NO.;
%let bad=1;
%end;
%if %upcase(&plot)^=T & %upcase(&plot)^=R & %upcase(&plot)^=K &
%upcase(&plot)^=N %then %do;
%put ERROR: PLOT MUST BE T, R, K, OR N.;
%let bad=1;
%end;
%let badalpha=0;
data _null_;
a=symget('alpha');
if a<0 | a>=1 then call symput('badalpha',1);
if a=0 then doalpha=0;
else doalpha=1;
call symput('doalpha',doalpha);
run;
%if &badalpha=1 %then %do;
%put ERROR: ALPHA MUST BE BETWEEN 0 AND 1;
%let bad=1;
%end;
%let nxs=0; %*count the number of x vars;
%do %until(%scan(&xvars,&nxs+1,' ')= );
%let nxs=%eval(&nxs+1);
%end;
%if &pvars^= %then %let plotvars=&pvars;
%else %let plotvars=&xvars;
%let nplotvar=0;
%do %until(%scan(&plotvars,&nplotvar+1,' ')= );
%let nplotvar=%eval(&nplotvar+1);
%end;
%do i=1 %to &nplotvar;
%let pv=%scan(&plotvars,&i,' ');
%let found=0;
%do j=1 %to &nxs;
%if &pv=%scan(&xvars,&j,' ') %then %let found=1;
%end;
%if &found=0 %then %do;
%put ERROR: PLOT VARIABLE &PV NOT ON XVARS LIST.;
%let bad=1;
%end;
%end;
%macro xlist(prefix); %*set-up var list code;
%do j=1 %to &nxs;
&prefix&j
%end;
%mend xlist;
%if &bad=0 %then %do;
data _setup; set &data; %*delete obs with mising values;
xmiss=0;
%do i=1 %to &nxs;
xx=%scan(&xvars,&i);
if xx=. then xmiss=1;
%end;
if &event=. or &time=. or xmiss=1 then delete;
%* run phreg and grab output datasets;
proc phreg data=_setup covout outest=_est ;
model &time*&event(0)= &xvars ;
%if &strata ^= %then %do;
strata &strata;
%end;
output out=_sch xbeta=xbeta ressch=rsch1-rsch&nxs
wtressch=wrsch1-wrsch&nxs;
data _sch1; set _sch(keep=&strata &time &event rsch1-rsch&nxs);
drop &event;
if &event=1;
rename
%do i=1 %to &nxs;
%let thisx=%scan(&xvars,&i,' ');
rsch&i=&thisx
%end;
;
proc sort; by &time;
data &outsch; set _sch1;
proc sort; by &strata &time;
** calculate overall Kaplan-Meier;
proc lifetest noprint data=_setup outs=_km;
time &time*&event(0);
data _km2; set _km;
keep &time probevt;
if _censor_=0; **keep event times;
probevt=1-survival;
label probevt='1-Overall Kaplan-Meier';
data _null_;
set _est;
%if %substr(&sysver,1,1)=6 %then %do;
if _type_='PARMS' & _name_='ESTIMATE';
%end;
%if %substr(&sysver,1,1)=8 %then %do;
if _type_='PARMS' & upcase(_name_)=upcase("&time") ;
%end;
%do i=1 %to &nxs;
%let thisx=%scan(&xvars,&i,' ');
call symput("est&i",&thisx);
%end;
data _sch2;
set _sch(keep=&strata &time &event wrsch1-wrsch&nxs);
keep &strata &time &xvars;
if &event=1;
%do i=1 %to &nxs;
%let thisx=%scan(&xvars,&i,' ');
&thisx=&&est&i + wrsch&i;
%end;
proc sort; by &time;
data _bt;
merge _sch2(in=ins) _km2(keep=&time probevt);
by &time;
if ins;
proc sort; by &strata &time;
run;
symbol1 i=join l=1 c=black;
symbol2 i=join l=2 c=black;
symbol3 i=join l=2 c=black;
symbol4 i=none v=dot h=.1 cm c=black;
%macro dovars(tvar);
%if %length(&tvar)=8 %then %let t=%substr(&tvar,1,7);
%else %let t=&tvar;
%daspline(&tvar,nk=&df,data=&outbt)
%do i=1 %to &nplotvar;
%let x=%scan(&plotvars,&i,' ');
data __temp1;
set &outbt;
&&_&t
%global _print_;
%let _print_=off;
%let ndummies=%eval(&df-2);
%if %substr(&sysver,1,1)=8 %then %do;
ods listing close;
%end;
proc genmod;
model &x=&tvar &t.1-&t.&ndummies/dist=normal obstats
%if &doalpha=1 %then %do;
alpha=&alpha
%end;
;
make 'OBSTATS' out=__temp2;
run;
%if %substr(&sysver,1,1)=8 %then %do;
ods listing;
%end;
data __temp3;
merge __temp1 __temp2;
/* lower=xbeta-2*std;
upper=xbeta+2*std; */
%if %upcase(&rug)=YES %then %do;
data __anno;
set __temp3;
x=&tvar;
y=0;
xsys='2';
ysys='1';
function='move';
output;
y=5;
function='draw';
output;
%end;
proc sort data=__temp3; by &tvar;
proc gplot data=__temp3 gout=gschbt
%if %upcase(&rug)=YES %then %do;
annotate=__anno
%end;
;
%if &doalpha=1 %then %do;
plot xbeta*&tvar=1
lower*&tvar=2
upper*&tvar=3
%if %upcase(&points)=YES %then %do;
&x*&tvar=4
%end;
/overlay vaxis=axis1 haxis=axis2
%if %upcase(&vref)=YES %then %do;
vref=0 lvref=3
%end;
;
%end;
%if &doalpha=0 %then %do;
plot xbeta*&tvar=1
%if %upcase(&points)=YES %then %do;
&x*&tvar=4
%end;
/overlay vaxis=axis1 haxis=axis2
%if %upcase(&vref)=YES %then %do;
vref=0 lvref=3
%end;
;
%end;
axis1 label=(r=0 a=90 "smooth(&x) df=&df");
%if &tvar=&time %then %do;
axis2 label=("&tvar");
%end;
%if &tvar=rtime %then %do;
axis2 label=("Rank for Variable &time");
%end;
%if &tvar=probevt %then %do;
axis2 label=('1-Overall Kaplan-Meier');
%end;
run; quit;
%end;
%mend dovars;
title2'Scaled residuals(Bt) as a fcn of time.';
proc rank data=_bt out=&outbt; var &time; ranks rtime;
proc corr pearson data=&outbt;
with &xvars; var &time rtime probevt;
label probevt='1 - Overall Kaplan-Meier';
%if %upcase(&plot)=T %then %do;
%dovars(&time)
%end;
%if %upcase(&plot)=R %then %do;
%dovars(rtime)
%end;
%if %upcase(&plot)=K %then %do;
%dovars(probevt)
%end;
run;
quit;
footnote; symbol; title2;
proc datasets nolist;
delete _setup _est _sch _sch1 _sch2 _km _km2 _bt
%if %upcase(&rug)=YES %then %do;
__anno
%end;
__temp1 __temp2 __temp3;
run;
quit;
%end;
%mend schoen;
